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I.  INTRODUCTION 


The  Advanced  Refractive  Effects  Prediction  System  (AREPS)  was  initially  developed 
as  a  performance  assessment  and  system  analysis  tool  for  U.S.  Navy  radar  systems 
operating  in  the  100  MHz  to  20  GHz  frequency  range.  Recently,  under  major  funding 
from  the  Office  of  Naval  Research  (ONR),  AREPS  capabilities  have  been  extended  to 
include  high  frequency  (HE,  2-30  MHz)  point-to-point  and  coverage  predictions  in 
support  of  military  communication  systems  and  a  fully  three-dimensional  ionospheric  ray 
tracing  capability  for  research  oriented  and  system  design  applications  at  HE.  Continuing 
this  effort  at  expanding  the  frequency  and  application  range  of  AREPS,  we  have  recently 
developed  a  medium  prediction  capability  for  coverage  analysis  in  the  AM  (300  kHz-  2 
MHz)  radio  band. 

With  these  additions  we  have  extended  the  AREPS  usable  prediction  frequency  range 
from  300  kHz  to  57  GHz  and  have  added  significant  communication  system  analysis 
capability.  In  this  document,  we  will  describe  the  latest  efforts  to  extend  the  range  of 
applications  of  AREPS  to  earth-satellite  communications.  The  new  capability,  to  be 
called  the  Earth-to-Satellite  Propagation  with  Meteorology  Model  (ESPM^),  employs  a 
ray  optics  propagation  model,  including  available  meteorological/refractivity  data,  to 
predict  signal  strength  at  satellite  locations.  Propagation  based  loss  mechanisms 
considered  in  the  signal  strength  determination  include  gaseous  absorption  and  absorption 
due  to  rain  on  the  propagation  path.  Since  we  intend  the  module  to  be  usable  for  both 
geo-stationary  and  orbiting  satellites,  we  have  also  included  an  orbital  prediction 
capability  which  allows  locations  and  time-of-availability  to  be  determined  for  almost 
any  satellite. 

This  document  is  intended  as  a  description  of  the  ray  optics  methods  which  have  been 
developed  for  use  in  signal  strength  prediction  in  ESPM^.  We  also  include  a  short 
discussion  (including  references  to  more  detailed  descriptions)  of  the  pre-existing  models 
chosen  for  prediction  of  losses  due  to  gaseous  and  rain  absorption,  as  they  directly  affect 
the  signal  strength  calculation,  for  completeness,  we  also  provide  a  short  discussion  of 
the  orbital  prediction  model  (with  references)  and  illustrate  its  usefulness  in  ESPM^.  In 
the  following  section  we  discuss,  generally,  the  refractive  model  that  is  used  in  the  ray 
optics  program.  Details  of  how  and  where  these  refractive  profiles  are  obtained  within 
AREPS  are  described  elsewhere. 

II.  ATMOSPHERIC  REFRACTIVITY  AND  MODEL  INPUT 

An  electromagnetic  wave  propagating  in  the  earth’s  atmosphere  feels  the  effects  of 
bound  electrons  in  the  molecular  species  which  make  up  the  atmosphere.  Under  the 
action  of  the  electric  field  of  the  passing  wave,  a  time  varying  induced  field  is  created 
which  radiates  a  secondary  field  that,  when  combined  with  the  original  wave,  may  cause 
the  path  of  the  original  wave  to  change.  This  change  in  the  original  path  is  called 
refraction  and,  under  certain  conditions,  may  cause  large  deviations  of  the  intended 
direction  of  the  original  wave.  For  the  wave  frequencies  used  in  satellite  communications 
and  the  nominal  molecular  density  of  the  lower  atmosphere,  it  is  usually  the  case  that  the 
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effects  of  refraction  can  be  ignored  and  the  effect  of  the  atmosphere  on  system 
performance  is  restricted  to  losses  due  to  rain  and,  to  a  lesser  extent,  gaseous  absorption 
of  wave  energy.  However,  the  refractive  effects  of  the  atmosphere  can  become 
anomalously  large  under  certain  atmospheric  conditions  of  temperature,  pressure  and 
humidity  (ducting)  and  the  refractive  effects  then  dominate  all  other  propagation  effects, 
essentially  determining  if  communication  is  possible  at  all.  Also,  when  a  satellite  appears 
low  on  the  horizon,  the  path  of  the  wave  through  the  atmosphere  is  very  long  and  the 
cumulative  refractive  effects  can  be  large  even  under  normal  refractive  conditions.  It  is 
for  these  anomalous  conditions  that  the  ESPM^  module  is  expected  to  be  most  useful. 

For  the  purposes  of  this  document,  we  assume  the  refractive  structure  of  the 
atmosphere  on  the  propagation  path  between  transmitting  antenna  and  receiving  antenna 
(assumed  throughout  this  document  to  be  at  the  satellite)  is  known.  In  general,  the 
atmosphere  varies  in  both  height  and  range  and  the  ray  optics  model  described  here 
allows  for  both.  To  include  the  range  dependence  in  a  consistent  way,  we  make  the 
following  definitions  regarding  the  propagation  path. 

Consider  the  plane  containing  the  transmitter  antenna  site,  the  instantaneous  receiver 
site  and  the  earth’s  center.  We  define  the  propagation  region  as  the  interior  portion  of  the 
plane  bounded  by  radials  from  the  earth’s  center  through  the  transmitter  and  receiver 
sites.  Within  this  region,  we  further  define  sectors  as  sub-regions  of  the  propagation 
region  within  which  the  earth’s  refractivity  can  be  considered  to  be  dependent  on  height 
only.  Formally,  let  ru  be  that  radial  within  a  sector  ‘k’  such  that  the  refractivity  changes 
most  ‘quickly’  as  we  move  from  the  sector  boundary  along  the  0-coordinate  within  the 
sector,  keeping  fixed.  For  this  fastest  changing  radial,  we  define  the  angle  \\ik  to  be  the 
angle  subtended  at  the  earth’s  center  over  which  the  refractivity  is  approximately 
constant.  The  angle  \^k  defines  the  angular  range  of  the  kth-sector.  Repeating  this  process 
throughout  the  propagation  region  for  sequential  sectors,  the  angular  range  of  the  region, 
i(/i,  which  is  also  the  angle  subtended  between  the  transmit  antenna  and  the  satellite,  can 
be  written 


K 

¥s  = 


k=\ 


(1) 


where  K  is  the  total  number  of  sectors  that  make  up  the  region.  For  each  sector  the 
refractivity  is  assumed  to  be  dependent  on  height  only.  We  note  that  this  formal  process 
for  defining  the  sectors  is  impractical  in  practice,  and  the  sectors  are  generally  defined 
based  on  the  availability  of  data  or  the  user’s  tolerance  for  model  execution  time. 

The  input  to  the  ray  trace  program  then,  consists  of  a  height  profile  of  refractivity 
(A(/),Z(i)),  i=  1,..,L,  where  Z(1)=0,  Z(L)=F(s  is  the  height  in  kilometers  of  the  satellite 
above  mean  sea  level  and  L  is  the  number  of  refractivity  samples,  for  each  of  the  sectors 
which  make  up  the  propagation  region.  For  heights  above  50  km  we  assume  the  earth’s 
atmosphere  can  be  well  approximated  as  vacuum  and  take  N{i)=Q  for  Z(/)  >  50  km.  For 
heights  below  50  km  we  further  assume  that  Z{i)  can  be  adequately  approximated  by  r{i), 
where  r  is  the  radial  coordinate  in  a  spherical  coordinate  system  referenced  to  the  earth’s 
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center.  So  at  all  heights  we  assume  we  can  make  the  approximation  that  ~  N(r(f)) 

where  r(i)=re+Z(z)  and  re  is  the  earth’s  radius  (K.D.  Anderson,  2000;  A.R.  Lowry,  et  ah, 

2002). 

We  also  require  specification  of  the  transmitter  height,  T,  in  kilometers  and  the  angle 
between  the  sub-transmitter  and  sub-satellite  points,  subtended  at  the  earth’s  center. 
The  angle  y/s,  which  can  be  alternately  specified  as  ground  range  between  the  sub¬ 
terminal  points,  defines  the  propagation  region.  The  height  of  the  satellite  Hs  is  also 
assumed  to  be  known  at  all  times  and  may  be  determined  by  the  satellite  orbit  prediction 
software  or  input,  if  known  to  acceptable  accuracy  by  the  user,  i.e.  geo- stationary 
satellite. 

The  calculation  of  signal  strength  requires  estimation  of  propagation  losses.  The 
model  included  for  estimation  of  the  gaseous  absorption  requires  input  of  the  height 
profile  of  the  absorption  rate.  In  practice,  this  profile  is  included  at  each  height  along  with 
the  refractivity.  The  rain  attenuation  loss  is  dependent  on  the  geographical  location  of  the 
transmit  antenna  and  so  we  also  require  the  latitude  and  longitude  coordinates  as  input. 

III.  RAY  CALCULATIONS 

I.  Basic  Ray  Equation  for  Range  Independent  Refractivity 

In  the  description  of  the  ray  tracing  program  within  this  section  we  will  assume  that 
the  transmitter  is  always  located  at  a  lower  height  than  the  receiver,  which  is  assumed  to 
be  the  satellite.  Given  the  symmetry  of  the  ray  tracing  method  we  could  interchange 
terminology  and  the  determined  rays  would  remain  unchanged.  However,  for  purposes  of 
this  discussion  we  will  use  ‘transmitter’  to  refer  to  the  lower  terminal  and  ‘receiver’  to 
refer  to  the  satellite. 

The  geometry  of  the  ray  trace  problem  for  a  propagation  region  consisting  of  a  single 
refractive  sector  (or  for  a  single  sector  in  a  multiple  sector  region)  is  shown  in  Figure  1. 
Since  the  refractivity  is  a  function  of  height  z  (or  radial  length  r)  only  within  a  sector,  the 
ray  will  remain  in  the  initial  azimuthal  plane,  which  we  take  to  be  the  (p  =0  plane.  In  Fig.. 
1,  an  initially  up-going  ray  is  launched  at  local  zenith  angle  6o  (<  90°)  from  a 
transmitting  antenna  located  at  the  point  (rT,0,0)  in  a  spherical  coordinate  system 
referenced  to  the  earth’s  center.  Here  rT=re+T  and  is  the  mean  earth’s  radius.  The 
problem  then  is  to  find  the  path  of  this  ray  through  the  specified  height-dependent 
refractivity  structure. 

In  Fig.  1  we  focus  on  the  ray  as  it  passes  through  a  spherical  atmospheric  layer 
bounded  by  r=r(i)=ri  and  r=r(i+l)=ri+\.  We  assume  the  ray  enters  the  layer  at  local 
zenith  angle  6i  and  passes  through  the  layer.  Then,  if  ds  is  an  infinitesimal  section  of  the 
ray  at  radial  distance  r  (n  <  r  <  ri+i)  making  a  zenith  angle  6  to  the  local  radial  (see  figure 
1),  we  can  write  (J.M.  Kelso,  1964) 


3 


Figure  1.  Ray  launched  at  angle  do  which  passes  through  the  layer  at  ri .  The  angle 
\|/s  is  the  angle  subtended  between  the  transmitter  at  (rr,0,0)  and  the 
satellite  located  at  (r5,\|/5,0). 


r 


dy/ 

dr 


=  tan  ^  = 


sin^ 

Vl  -sin^  9 


(2) 


where  dy/  is  the  angle  subtended  at  the  center  of  the  earth  by  ds. 

We  next  apply  Bouger’s  Law  (Snell’s  Law  equivalent  for  spherical  geometry),  which 
states  that 


n{r)r  sin  9  =  n{r.  )r.  sin  0^  =  n^rj  sin  6^ , 


(3) 


in  (2)  to  obtain 


_  u,?;.sin^, 

sin^  0. 


(4) 


Here  n(r)  is  the  index-of-refraction  at  the  radial  distance  r  which  is  obtained  from  the 
refractivity  N(r)  by  the  relation 


Nir)  =  inir)- 1)10\ 


(5) 
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and  ni=n(ri)  is  the  index  at  the  bottom  of  the  layer.  We  then  obtain  the  total  angle 
subtended  by  the  ray  in  passing  through  the  layer  by  integrating  (4) 


n.sind. 
r.  ^  r 

I  r: 


n<r<  r\+\  (6) 


Equation  (6)  is  the  basic  ray  equation  evaluated  in  the  ray  trace  model. 

2.  Evaluation  of  the  Basic  Ray  Equation 

One  of  our  goals  in  solving  the  basic  ray  equation  (6)  is  minimization  of  the 
computation  time  required  to  obtain  a  solution.  This  is  especially  important  in  light  of  the 
homing  procedure  to  be  discussed  below  which  generally  requires  many  evaluations  of 
the  basic  ray  equation  to  obtain  a  solution  for  a  small  error  tolerance.  Minimization  of 
computation  time  can  be  partially  accomplished  by  making  simplifying  assumptions  in 
the  integrand  in  order  to  cast  (6)  into  a  form  for  which  closed  form  evaluation  is  possible. 

The  first  simplification  of  (6)  involves  an  approximation  for  the  term  (  rjr)  which 
appears  twice  in  the  equation.  We  write  this  term  as, 


r:^2 


{-y  =  ( 

r  r+z 


^-y  =  (- 


+  z, 


r  +  z,  +  z  -  z, 


-)^  = 


{i+hf 


where  h=(z-zi)/ri  and  Zi  <  z  <  Zi+i.  We  assume  that  h  is  much  smaller  than  unity  at  all 
atmospheric  heights  and  expand  this  term  to  obtain  the  approximation 

{^y  ^l-2h  +  3h^ - .  (7) 

r 

We  further  simplify  the  ray  equation  by  assuming  that  the  refractivity  is  a  linear 
function  of  height  in  the  region  between  the  layer  boundaries  and  write  N{r)=Ni  +  S{{r-r^, 
where  n  <  r  <  n+i  and  Si  is  the  constant  refractivity  slope  in  the  i-th  layer  given  by 
Si={Ni+\-N^I{ri+\-r^.  From  (5)  we  can  write  the  index-of-refraction  in  the  layer  as 

n{r)  =  n.  +5,.(r-r.) 

'  '  '  (8) 

=io-®7v.  +1,.?.  =io-®«s. 

so,  since  (6)  actually  requires  the  square  of  the  index-of-refraction,  we  have 
n  {h)=  Hi  +  Isifiirih  +  Si  ri  h  where  and  rii  =2ni-\. 

With  these  simplifications  the  basic  ray  equation  (6)  can  be  written  as 
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¥\^:'=c, 


2  3 

r.— r,+-r3 

r  r. 


where  Co=«isin6*i  and  we  have  the  following  definitions, 

''p  hdh  ^ 

J  /  2 

Q^a  +  bh  +  ch 

where  hy^{z{+\-z\)lri,  a=ni  cos  ,  b=2ni5iri  +  2ni  sin  Oi  and  c=Si  n  -  3ni  sin  Oi.  The 
integrals  T2  and  Ts  can  be  written  in  terms  of  Ti^  which  is  a  standard  form  found  in 
integral  tables  (R.S.  Burington,  1965).  Then,  after  some  algebraic  manipulation  (see 
Appendix),  we  obtain  the  following  closed  form  expression  for  the  angle  subtended  by 
the  ray  in  passing  through  the  atmospheric  layer  at 

w  1;™  =  Q  \a(F{IiJ-  Fm  +  (^a  +  bh,+clilBOi, )  -  ^/^^(0))J  (10) 

where  for  c  <  0, 


h^dh 


’J. 


a  +  bh  +  ch^ 


(9) 


r.  =  J 


dh 


■\Ja  +  bh  +  ch 


.,r,= 


F{h)  = 


sin  '(a(h)) 


,a{h)  = 


-  2ch  -  b 
ylb^-4ab  ’ 


and  for  c  >  0, 


F(h)  =  -^log(2ch  +  b  +  2-\[c^ja  +  bh  +  ch^  ) . 

Vc 

Here  B{h)=3h/2cr^  -  2/cri  -  9b/4c^ri^  and  ^4=1  +  b/cri  +  9b^/8c^ri^  -  3aJ2cn  ■  The  result 
for  the  special  case  c=0,  along  with  the  modifications  necessary  for  a  down-going  ray,  are 
given  in  the  Appendix. 

The  expression  (10)  is  evaluated  for  each  atmospheric  layer  and  the  total  angle 
subtended  by  the  ray  within  the  sector,  \\fk,  is  obtained  from 


(=1 

where  T*  is  the  subset  of  height  layers  traversed  by  the  ray  within  sector  k. 

In  the  homing  process,  the  sum  of  the  \|/^‘s  is  calculated  for  a  given  ray  and  is 
compared  to  the  angle  subtended  by  the  sub-satellite  point,  ii/^.  If  the  subtended  angles 
differ  by  less  than  some  preset  tolerance  level  we  take  that  ray  as  the  homed  ray  and  the 
initial  launch  angle  for  the  ray,  9o,  is  saved.  If  the  difference  between  subtended  angles 
exceeds  the  tolerance  level  a  new  ray  with  a  different  initial  angle  is  launched  and  the 
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process  continues  until  the  desired  ray  is  found.  We  discuss  the  homing  procedure  in 
more  detail  in  a  later  section.  In  the  next  section  we  consider  scenarios  in  which  rays  may 
be  reflected  in  atmospheric  layers. 

3.  Ray  Turning  Points  in  the  Atmosphere  (Ducting) 

In  the  previous  sections  we  described  the  general  method  for  determining  rays  in  the 
atmosphere.  It  was  implicitly  assumed  throughout  that  the  rays  were  launched  from  the 
transmitting  antenna  at  initial  zenith  angles  6q  <  90°  and  that  the  vertical  gradients  of 
refractivity  were  weak  in  the  sense  that  such  initially  ‘up-going’  rays  remained  up-going  , 
i.e.  the  local  zenith  angles  6*i  <  90°  for  all  /=!,. . .,  L.  In  general,  however,  an  initially 
up(down)-going  ray  may  reverse  direction  at  ray  ‘turning  points’  in  the  atmosphere.  In 
this  section  we  discuss  these  turning  points  and  describe  the  method  for  dealing  with 
them  in  the  ray  trace  program. 

a.  Preliminaries 


We  begin  by  examining  the  effect  of  vertical  refractive  gradient  on  ray  propagation 
within  the  atmosphere.  Consider  again  the  scenario  shown  in  Fig.  1  and  focus  on  the 
refractive  layer  between  n  and  ri+\.  From  equation  (3)  we  can  write 


sin  =  fi.  sin  0. 
n  r. 

A  -(—)(—) 


(11) 


Since  n+i  >  ri  the  layer  constant  Pi  will,  in  general,  be  less  than  or  greater  than  unity 
depending  on  the  change  of  the  index  of  refraction  in  the  layer.  In  the  earth’s  lower 
atmosphere,  the  index  of  refraction  is  generally  a  decreasing  function  of  height  so  we 
have  Ui+i  <  where  equality  holds  for  a  uniform,  isotropic  atmosphere  and  by  a  ‘weak’ 
refractive  gradient  we  mean  that  ni/ui+i  ~  1. 

Consider  first  the  scenario  discussed  earlier  which  we  have  termed  the  ‘weak’  vertical 
refractive  gradient  situation.  In  this  case  the  constant  Pi  is  less  than  unity  for  all  layers 
L  and,  from  equation  (1 1),  we  see  that  the  local  zenith  angle  at  the  top  of  a  layer  is 
always  less  than  that  at  the  bottom  of  the  layer.  Thus  a  ray  that  starts  out  up-going  (^i  < 
90°)  will  remain  up-going  as  it  propagates  through  the  atmosphere  and  we  can  define  a 
weak  vertical  refractive  gradient  as  <  1  or  njni+i  <  ri+iM. 

Next  consider  a  ray  which  is  initially  down-going  as  it  propagates  through  the  weakly 
refracting  atmosphere.  Maintaining  our  height  ordering  such  that  n+i  >  ri,  we  have 


sin^,.  =  sin^.^i 
=  {—){—) 


(12) 
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where  6i+\  >  90°.  Pi+i  is  again  greater  or  less  than  unity  depending  on  the  refractive 
gradient  in  the  layer.  In  this  case,  for  a  weak  refractive  gradient,  Pi+\  >  1  for  all  layers, 
the  local  zenith  angle  on  exiting  the  layer,  6{,  is  smaller  than  the  incident  angle  and  we 
again  have  «i+i/«i  >  r-Jn+x,  However,  if  the  local  incident  zenith  angle  is  very  near  90°  for 
some  layer,  the  sine  of  the  exit  angle  may  exceed  unity  and  become  undefined.  In  this 
case  a  ray  turning  point  has  been  reached  within  the  layer  and  the  ray  trace  is  stopped. 
The  local  zenith  angle,  now  less  than  90°,  is  determined  as  described  below  and  the  ray 
trace  is  re-started  from  that  point.  These  scenarios  for  weak  refractive  gradients,  for  both 
up-going  and  down-going  initial  rays,  are  shown  in  the  top  panel  of  Figure  2. 

Next  consider  a  scenario  in  which  at  least  one  atmospheric  layer  has  a  strong 
refractive  gradient.  Then,  from  equation  (1 1)  for  an  up-going  ray,  Px  >  1  for  such  a  layer 
and  the  local  zenith  angle  of  a  ray  exiting  at  the  top  of  the  layer  is  greater  than  that  at  the 
bottom  of  the  layer.  Additionally,  if  the  local  zenith  angle  Ox  at  the  bottom  of  the  layer  is 
less  than  but  close  to  90°,  the  sine  of  the  exiting  angle  may  exceed  unity  and  become 
undefined.  In  this  instance  the  ray  has  turned  within  the  layer  and  again  the  ray  trace 
program  must  deal  with  the  change  in  the  ray  trajectory. 
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Figure  2.  Illustrating  the  relationship  between  layer  incident  and  exit  angles  for 
(top)  a  ‘weak’  refractive  layer,  and  (bottom)  for  a  ‘strong’  refractive 
layer.  Arrows  indicate  the  direction  the  zenith  angle  moves  as  the 
ray  passes  through  the  layer  for  each  situation. 

For  an  initially  down-going  ray  incident  on  a  layer  for  which  Pi+\  <  1,  the  local  zenith 
angle  on  exiting  the  layer  will  exceed  the  incident  angle.  If  the  layer  is  near  the  earth’s 
surface,  the  ray  may  hit  the  earth  and  the  ray  trace  is  stopped.  These  large  refractive 
gradient  scenarios  are  shown  in  the  bottom  panel  in  Fig.  2  for  both  up-going  and  down¬ 
going  initial  rays. 

b.  Ray  methods  at  turning  points 

In  the  previous  section  we  showed  that  a  ray  may  reverse  direction  in  a  refractive 
layer  under  some  conditions.  To  determine  if  such  a  ray  turning  occurs  in  a  layer  we 
apply  equation  (1 1)  for  initially  up-going  rays  or  equation  (12)  for  down-going  rays  to 
each  layer.  If,  for  an  up-going  ray  (^i  <  90°)  incident  on  the  layer,  the  local  zenith  angle  at 
exit  (^i+i)  exceeds  90°,  the  ray  has  turned  at  some  height  rgo  at  which  0=90°,  where  n  < 

^90  <  n+i-  We  can  determine  rgo  from  equation  (I  I)  by  noting  that  at  the  turning  height 
we  have 


1  =  (— x— )sin^,, 


^90  ^90 


(13) 
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where  n^Q=n{  +  Si{r^Q-r^.  Substituting  this  expression  for  in  equation  (13)  we  find  that 
rgo  satisfies  a  quadratic  equation  which  we  solve  for  the  turning  height.  We  then  restart 
the  ray  trace  program  at  that  height,  adjusting  parameters  for  the  now  oppositely  directed 
ray.  Similar  results  apply  for  the  initially  down-going  ray  which  turns  ‘up-ward’  in  the 
layer.  Figure  3  shows  a  vertical  fan  of  rays  launched  into  a  strong  refractive  duct  (shown 
in  the  figure).  As  indicated  in  Fig.  2,  the  ray  behavior  in  the  duct  depends  critically  on 
both  the  incident  ray  angle  and  the  refractive  gradient  in  the  layer. 


4.  Homing 

We  indicated  earlier  that  the  determination  of  the  ray(s)  which  connects  transmitter 
and  satellite  requires  a  homing  procedure.  If  the  satellite  lies  above  the  geometric  horizon 
of  the  transmitter  there  will,  in  general,  be  two  rays  which  connect  transmitter  to  receiver, 
the  direct  ray  and  the  earth  reflected  ray.  However,  the  existence  of  the  reflected  ray  also 
requires  that  the  transmitting  antenna  be  elevated  above  the  ground.  In  this  section  we 
describe  the  homing  procedure  implemented  in  this  version  of  the  program. 

a.  Direct  ray 

By  direct  ray  we  mean  that  ray  which  propagates  directly  from  transmitter  to  receiver 
with  no  intervening  earth  or  structural  interaction.  This  mode  will  generally  exist  if  the 
satellite  is  above  the  geometric  horizon  of  the  transmitting  antenna.  It  may  also  exist  for  a 
satellite  slightly  below  the  horizon  if  a  sufficiently  strong  refractive  gradient  exists  in  a 
non-uniform  atmosphere.  This  ‘beyond-line-of-sight’  mode  does  not  exist  for  a  range 
independent  propagation  region  since  a  trapped  ray  remains  trapped  if  there  is  no 
horizontal  gradient  in  the  refractivity. 
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RAY  FAN 


REFRACTIVITY  PROFILE 


Figure  3.  Ray  fan  launched  into  a  strongly  refracting  duct.  The  transmitter  lies 
within  the  duct  and  all  launch  angles  are  near  90°.  Ray  behavior  in  the 
duct  is  critically  dependent  on  launch  angle  and  the  strength  of  the  duct. 

The  homing  procedure  for  the  direct  ray  proceeds  by  a  straight-forward  iterative 
procedure  involving  the  launching  of  rays  at  prescribed  launch  angles  to  the  height  of  the 
satellite,  using  the  ray  equations  described  above,  and  comparing  the  resultant  subtended 
angles  (or  equivalently,  the  ground  range  between  the  sub-transmitter  and  sub-ray 
terminal  points)  to  the  angle  (ground  range)  between  the  transmitter  and  sub-satellite 
point.  The  direct  ray  is  taken  to  be  that  ray  for  which  the  difference  is  less  than  some 
prescribed  tolerance  parameter.  This  procedure  has  the  important  advantage  of 
robustness,  in  that  a  solution  can  almost  always  be  found. 

The  execution  time  required  to  locate  the  direct  ray  is  determined  by  the  number  of 
iterations  required  to  satisfy  the  specified  tolerance  which,  for  the  current  version  is  .01 
meters  in  ground  range.  Clearly,  a  good  initial  estimate  of  the  required  ray  launch  angle 
is  crucial  in  reducing  this  execution  time.  Generally  speaking,  for  most  transmission 
angles  and  refractive  conditions,  the  direct  ray  differs  very  little  from  the  straight  line 
path  between  the  transmitter  and  satellite.  For  this  reason  the  initial  search  angle  used  in 
the  model  is  the  angle  corresponding  to  this  straight  line  path. 
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The  procedure  begins  by  launching  a  ray  at  the  straight  line  zenith  angle  and  tracing 
to  the  satellite  height.  If  the  difference  between  the  required  ground  range  and  the 
spanned  ground  range  is  greater  than  zero  the  initial  launch  angle  is  increased  by  an 
initial  increment  d9=2)°  and  a  new  ray  is  launched.  If  the  initial  ground  range  difference  is 
less  than  zero  the  initial  launch  angle  is  decreased  by  dO  for  the  next  ray.  Each  time  the 
sign  of  the  range  difference  changes  the  launch  angle  increment  is  halved  and  the  sign  of 
the  launch  angle  increment  is  reversed.  This  procedure  continues  until  the  ground  range 
difference  falls  below  the  tolerance  level. 

An  illustration  of  a  typical  homing  procedure  for  the  direct  ray  is  shown  in  Figure  4. 
We  note  the  relatively  large  swings  in  the  difference  between  the  required  ground  range 
and  the  candidate  ray  ground  range.  These  initially  large  swings  result  from  the  fact  that 
the  procedure  requires  two  estimates  of  range  difference  to  establish  the  next  launch 
angle.  The  initial  estimate  of  the  straight  line  launch  angle  is  incremented  by  a  fixed  +1° 
or  -1°  increment  for  the  second  iteration,  which  provides  the  estimate  for  the  third 
iteration.  This  fixed  1°  increment  does  not  necessarily  result  in  a  ‘good’  estimate,  as 
indicated  by  the  scenario  shown  in  Fig.  4. 

b.  Earth  reflected  ray 

Homing  for  the  reflected  ray  is  significantly  more  complicated  in  that  it  requires  a 
two-dimensional  iteration  scheme  since  the  location  of  the  reflection  point  on  the  earth’s 
surface  is  not  known  initially.  The  requirement  for  the  location  of  the  reflection  point  is 
that  the  reflection  angle  of  the  ray,  measured  relative  to  the  plane  tangent  to  the  earth  at 
the  reflection  point,  is  equal  to  the  incident  angle.  This  is  an  approximation  to  the  full 
theory  of  the  reflection  process  from  a  curved  surface  (reference)  which  becomes 
increasingly  inaccurate  for  very  low  incident  angles.  We  use  the  approximation  here 
because  it  provides  a  fast,  acceptably  accurate  estimate  for  the  homing  process  in  most 
cases. 

The  process  begins  by  selection  of  an  initial  candidate  point  for  the  reflection  location 
on  the  earth’s  surface.  In  most  cases,  the  actual  location  of  the  reflection  point  is  close  to 
the  transmitter  and  will  generally  lie  within  the  refractivity  sector  nearest  to  the 
transmitter  site.  For  the  purposes  of  the  homing  procedure,  the  selected  candidate 
reflection  point  is  treated  as  the  transmitter  location  and  direct  rays  connecting  this  point 
to  the  satellite  and  actual  transmitter  location  are  determined  by  the  method  described  in 
the  previous  section.  The  launch  angles  of  the  determined  rays  are  compared  and  the 
candidate  location  of  the  reflection  point  is  adjusted  based  on  the  difference  in  the  angles. 
This  process  continues  until  a  point  is  determined  for  which  the  difference  between 
launch  angles  for  the  direct  rays  to  the  actual  transmitter  site  and  satellite  is  less  than  a 
specified  tolerance. 
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Satellite  Height=1000  km 
Transmitter  Height=100  m 
Satellite  Offset=20°(SSd=2226.39  km) 
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Figure  4.  Illustration  of  the  homing  process  for  a  typical  direct  ray  determination. 
Each  point  represents  a  candidate  ray  launched  at  a  different  initial 
angle.  The  vertical  axis  is  the  difference  between  SSd  (=  tf/s=20°),  the 

fixed  ground  distance  between  the  sub-satellite  and  sub-transmitter 
points,  and  the  ground  range  subtended  by  the  candidate  ray  at  each 
iteration,  Rd. 


The  result  of  this  homing  process  is  the  identification  of  that  ray  which  connects 
transmitter  to  satellite.  For  later  processing  we  also  require  the  total  ray  length,  x,  of  this 
ray.  The  total  path  length  is  obtained  by  summing  the  individual  path  length  increments, 
Xi,  from  each  refractive  layer.  From  Fig.  1  we  see  that  the  path  length  increment  in  the 
ith-  layer  is  given  by 


cos(^^  ) .  (14) 

5.  Modifications  for  Range  Dependent  Ray  Calculations 

It  has  been  assumed  in  all  of  the  above  that  the  refractivity  varies  in  height  only.  In 
reality,  the  earth’s  refractivity  also  varies  in  the  horizontal  direction  (range)  and  some 
knowledge  of  that  variation  is  required  for  accurate  determination  of  rays,  especially  for 
the  ducting  conditions  described  earlier.  If  we  consider  the  instantaneous  vertical  plane 
which  contains  the  transmitter  location,  satellite  location  and  the  earth’s  center,  we  earlier 
defined  the  propagation  region  as  that  portion  of  the  plane  bounded  by  radials  from  the 
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earth’s  center  through  each  of  the  satellite  and  transmitter  locations.  Within  the 
propagation  region,  it  is  generally  true  that  the  earth’s  refractivity  at  a  fixed  height 
changes  as  we  move  within  the  region.  We  defined  a  sector  as  that  portion  of  the  region 
for  which  the  refractivity  is  approximately  constant  as  we  move  in  range.  Thus  a 
propagation  region  may  be  made  up  of  several  sectors,  within  each  of  which  the 
refractivity  varies,  approximately,  in  height  only  and  so  all  of  our  previous  results  are 
assumed  to  be  valid  within  a  sector.  If  a  ray  crosses  into  different  sectors  as  it  progresses 
to  the  satellite,  we  must  modify  the  ray  trace  to  reflect  the  different  refractive 
environment  in  each  new  sector.  In  this  section  we  describe  this  process. 

Consider  the  scenario  shown  in  Figure  5,  where  an  up-going  ray  is  shown  incident  on 
a  layer  at  height  where  the  superscript  denotes  the  ‘k-th’  sector  and  is  the  total 
angle  subtended  by  the  sector.  Several  scenarios  are  possible  for  the  ray  within  the  layer. 
First,  if  the  ray  does  not  reach  a  turning  point  in  the  layer  and  the  total  angle  subtended  in 
the  layer  plus  the  sum  of  the  angles  subtended  in  previous  layers  does  not  exceed  we 
continue  the  ray  to  the  next  layer.  Second,  if  a  turning  point  is  reached  in  the  current  layer 
we  deal  with  the  reflection  as  described  earlier.  In  this  case  we  sum  the  angles  subtended 
before  and  after  reflection  within  the  layer  and  add  this  to  the  total  angle  subtended  in  the 
previous  layers.  Again,  if  the  total  does  not  exceed  \|/fo  the  ray  trace  is  continued  with  the 
ray  now  going  in  the  opposite  direction. 

Next  we  assume  that  the  ray  does  not  reach  a  turning  point  in  the  initial  layer  but  it 
does  pass  into  the  next  sector  at  the  transition  height  (see  Fig.  5).  To  determine  the 
transition  height  we  use  an  iterative  procedure  of  selecting  a  candidate  height  between 
and  /i+i  and  calculating  the  angle  subtended  by  the  ray  in  traversing  the  ‘virtual’  layer 
between  and  the  candidate  transition  height.  The  subtended  angle  in  the  virtual  layer  is 
added  to  the  total  angle  subtended  in  the  lower  layers  (up-going  initial  ray)  and  the  result 
is  again  compared  to  Based  on  the  comparison,  the  candidate  transition  height  is 
adjusted  up  or  down  and  the  process  repeated  until  the  difference  between  the  total  angle 
subtended  by  the  ray  up  to  the  transition  height  and  is  less  than  some  pre-set  criteria. 
Then,  applying  Bouger’s  Law,  the  local  zenith  angle  of  the  ray  at  rt  is  calculated  and  the 
ray  is  continued  at  that  angle  into  the  next  sector.  The  ray  trace  then  continues  as  before, 
using  the  refractivity  profile  assigned  to  the  following  sector. 

The  final  possibility  involves  a  situation  in  which  the  ray  not  only  passes  into  the  next 
sector  within  a  layer,  but  also  reaches  a  turning  point  within  the  layer.  In  this  case  we 
must  calculate  the  turning  point,  rpo,  and  transition  height,  Vt,  within  the  layer  to 
determine  which  height  was  reached  first.  If  the  turning  point  is  reached  before  the 
transition  height,  the  turning  occurs  before  the  ray  passes  into  the  next  sector  and  the 
turning  of  the  ray  is  handled  in  the  usual  way  within  the  current  layer/sector.  The  ray, 
now  propagating  in  the  opposite  direction  is  then  followed  into  the  next  sector  as 
described  in  the  previous  paragraph. 
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Figure  5.  Illustrating  ray  transition  between  adjacent  propagation  sectors.  The 
radial  of  length  rt  (=  Vg+hi)  is  the  transition  height  between  sectors.  The 
angles  \(/a  and  \(/a+/  are  the  angles  subtended  by  each  sector. 


If  the  transition  to  a  new  sector  occurs  at  a  lower  height  than  the  ray  turning  point,  the 
transition  height  is  determined  as  described  above  and  the  ray  continued  in  the  next 
sector.  The  ray  may  or  may  not  reach  a  turning  point  in  the  new  sector  where  a  different 
refractivity  profile  exists.  If  it  does  reflect  in  the  new  sector,  the  usual  methods  are 
employed  within  the  new  refractive  profile  for  the  next  sector. 

Figure  6  shows  a  scenario  where  the  propagation  region  consists  of  two  sectors.  The 
refractivity  in  the  initial  sector  contains  a  trapping  layer  while  the  second  does  not.  As 
indicated,  the  initial  fan  of  rays  are  all  trapped  in  the  duct  (top  right  of  Fig.  6).  As  they 
propagate  to  the  next  sector  the  rays  are  no  longer  bound  to  the  low  altitude  duct  and  may 
escape  to  reach  satellite  heights. 
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Figure  6.  Illustrating  the  effect  on  ray  propagation  of  range  dependent  refractivity. 
The  initial  sector  contains  a  ducting  profile  which  traps  rays  launched 
within  the  duct  (upper  right).  When  the  trapped  rays  pass  into  the  next 
non-ducting  sector  they  can  escape  to  satellite  heights  (lower  right). 


IV.  SATELLITE  LOCATION  AND  TRACKING 

Explicit  in  the  previous  discussion  regarding  homing  is  the  assumption  that  the 
location  (Rs,^s)  of  the  satellite  is  known  whenever  one  wishes  to  communicate.  While 
this  assumption  may  be  partially  true  when  dealing  with  geo-stationary  satellites,  it  is 
certainly  not  true  for  low  orbiting  satellites  whose  location  changes  rapidly  in  time.  Such 
satellites  may  be  available  for  transmissions  from  a  fixed  transmitter  site  for  only  a  few 
minutes  several  times  per  day  and  so  prior  knowledge  of  those  occurrences  can  be  critical 
for  successful  system  operation. 

The  accurate  prediction  of  the  location  of  a  satellite  at  any  given  time  requires,  first, 
precise  knowledge  of  its  location  and  orbital  velocity  at  some  reference  time  as  near  to 
the  required  time  as  possible.  This  precise  information  is  used  as  an  initial  condition  for 
the  solution  of  the  orbital  equation  of  motion,  which  solution  provides  the  time 
development  of  the  orbit  and  allows  future  (and  past)  satellite  locations  and  velocities  to 
be  determined.  In  time,  however,  inaccuracies  in  this  process  (caused  by  inaccurate  input 
and/or  incomplete  knowledge  of  the  forces  acting  on  the  satellite)  cause  the  actual 
satellite  location  and  velocity  to  deviate  from  the  predicted  values.  When  the  difference 
between  actual  and  predicted  satellite  location  and  velocity  becomes  excessive,  a  new 
precise  set  of  input  data  must  be  obtained  and  the  process  of  determining  the  time 
dependence  of  the  orbital  parameters  is  repeated. 
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To  determine  satellite  location  for  the  ESPM^  we  have  included  two  orbital  prediction 
models,  designated  SGP4  and  SDP4,  developed  for  the  North  American  Aerospace 
Defense  Command  (NORAD),  which  is  a  ‘bi-national  United  States  and  Canadian 
organization  charged  with  the  missions  of  aerospace  warning  and  aerospace  control  for 
North  America.  Aerospace  warning  includes  the  monitoring  of  man-made  objects  in 
space,  and  the  detection,  validation,  and  warning  of  attack  against  North  America 
whether  by  aircraft,  missiles,  or  space  vehicles,  through  mutual  support  arrangements 
with  other  commands.  Aerospace  control  includes  ensuring  air  sovereignty  and  air 
defense  of  the  airspace  of  Canada  and  the  United  States’  (http://www.norad.mil).  To 
accomplish  this  mission,  NORAD  maintains  and  publishes  the  precise  satellite  location 
information  necessary  for  the  determination  of  satellite  location  for  all  resident  space 
objects  (F.R.  Hoots  and  R.L.  Roehrich,  1988).  This  information  is  available  on-line  in  the 
form  of  two-line  element  sets  (TLE)  for  each  satellite.  The  TEE  provides  the  input  data 
for  execution  of  the  orbital  prediction  models. 

The  models  SGP4,  for  near  earth  satellites  (orbital  period  less  than  255  minutes)  and 
SDP4  for  deep  space  orbiters,  including  geo- stationary  satellites,  require  NORAD  TLE 
orbital  element  fdes  for  execution.  The  driver  program  selects  the  appropriate  model 
based  on  the  input  data,  a  process  that  is  transparent  to  the  user.  The  basic  output  from 
the  program  consists  of  the  satellite  location  and  velocity  at  the  required  time  in  the 
inertial  coordinate  system  used  for  the  calculations.  Documentation  of  the  details  of  the 
programs  SGP4  and  SDP4  are  available  (F.R.  Hoots  and  R.L.  Roehrich,  1988;  and 
references  therein)  and  are  not  presented  here  as  they  are  beyond  the  scope  of  this 
document. 

For  the  purposes  of  the  ESPM^  program,  software  has  been  developed  for  translation 
of  the  satellite  location  to  the  earth-centered  spherical  coordinate  system  which  rotates 
with  the  earth  that  is  used  for  ray  determination  (T.S.  Kelso,  1995(a);  1995(b);  1996). 
Once  the  location  is  determined  in  the  proper  coordinate  system,  we  can  calculate  the 
location  (i?s,Vs)  required  for  the  ray  trace  homing  procedure.  We  also  provide  prediction 
of  the  times  at  which  the  satellite  climbs  above  the  local  horizon,  and  is  thus  available  for 
communication  from  a  given  transmitter  site  for  a  non-ducting  refractive  profde.  Figure 
7  shows  example  output  which  combines  the  satellite  tracking  software  with  ESPM^  to 
provide  estimates  of  signal  strength  at  the  satellite  and  transmit  antenna  pointing  angles. 

It  should  be  stressed  that  a  user  with  independently  obtained  satellite  location  data  is  not 
required  to  use  the  satellite  orbit  programs  described  here.  However,  such  location  data 
must  be  entered  in  a  format  consistent  with  that  required  by  the  ESPM^  software. 
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GCA°=  Great-Circle  Angle  in  degrees 

T°=  Pointing  Angle  at  Transmitter  in  degrees  from  zenith 

dB(V/m)=  deci-Bels  relative  to  1  v/m  for  1  W  transmitted 


Figure  7.  Illustration  of  the  use  of  the  satellite  orbit  determination  software  with 
the  ESPM^  ray  trace  software.  Green  dots  represent  the  sub-satellite 
locations  for  an  orbital  pass  of  the  LANDSAT  5  satellite.  The  homed 
rays  are  shown  in  blue.  The  table  contains  estimates  of  signal  strength 
and  antenna  pointing  angles  for  rays  transmitted  from  the  assumed 
transmitter  site. 


V.  SIGNAL  STRENGTH  EVALUATION 

In  the  preceding  sections  we  described  the  ray  tracing  process  and  the  homing 
procedure  which  permits  the  ray(s)  which  connect  transmitter  and  receiver  to  be 
determined.  With  this  information  in  hand,  we  can  proceed  to  the  evaluation  of  the 
electric  field  strength  at  the  receiver  in  order  to  assess  the  expected  system  performance. 
In  this  section  we  derive  general  equations  for  the  fields  and  discuss  the  models  used  for 
evaluation  of  the  terms  which  appear  in  the  equations. 

I.  General  Field  Equations 

We  assume  that  the  total  vector  field  at  the  satellite  position  can  be  written, 

Et  =  Ed  )  +  Er  ) ,  (15) 
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where  |£'t|  is  the  total  field,  |£'d|  is  the  direct  field  amplitude,  \Er\  is  the  earth  reflected  field 
amplitude,  k  is  the  free  space  wave  number  and 

<t>r  =  +arg(i?)/A: 

Here  /)o  is  the  initial  phase  of  the  direct  ray  relative  to  the  reflected  ray  and  pd(r)  is  the 
phase  path  length  given  by. 


5  g  s 

=^nds^,p^  =^nds^+^nds^  ,  (17) 

'  t  g 

where  t  is  the  transmitter  location,  5  is  the  satellite  location  and  g  is  the  ground  reflection 
point.  In  this  equation,  dsd(r),  is  a  ray  element  along  the  direct  (reflected)  ray  path,  which 
is  determined  in  the  ray  trace  program,  and  n  is  the  index  of  refraction.  The  method  used 
to  determine  the  phase  path  is  entirely  similar  to  that  used  for  the  rays  and  is  outlined  in 
the  Appendix.  R  is  the  complex  Fresnel  reflection  coefficient  which  accounts  for  the 
effects  of  a  finitely  conducting  earth,  and  its  formula  is  provided  in  the  Appendix  for  both 
horizontal  and  vertical  polarization. 

The  reflected  field  amplitude  in  (15)  is  composed  of  several  terms  to  account  for 
reflections  from  a  non-perfectly  conducting,  curved-earth  surface.  Thus  we  assume  the 
reflected  field  amplitude  can  be  written. 


Er=Eo\R\D,  (18) 

where  |£’o|  is  the  amplitude  of  a  wave  reflected  from  a  smooth,  flat,  perfectly  conducting 
tangent  plane  at  the  reflection  point  after  traversing  a  total  distance  Jr  to  the  satellite.  D  is 
the  (real)  divergence  of  the  wave  due  to  reflection  from  a  curved  earth  surface. 

From  (15)  the  magnitude  of  the  total  field  is  given  by 

\Et[  =  Qxp{ik(l)^ )  +  E^  exp(/A:^^, )){E^  exp(/A:<zi^ )  +  E^  exp(/A:^^, )) * , 

where  the  asterisk  represents  complex  conjugate.  We  set  /io=0  and  assume  all  the  fields 
are  co-polarized.  Then  the  field  amplitudes  can  be  taken  as  real  scalars  and  the  total  field 
at  the  satellite  is  given  by 


=  E]  +E;  +2E,E,  cosW,  (19) 


where  the  phase  of  the  reflected  wave  includes  the  Fresnel  phase  as  indicated  in  (16). 
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Now  the  power  flux  density  at  the  satellite  in  the  direct  ray,  IZq^  is  given  from  the 
Poynting  Theorem  by 


^  =  Pt.gJsKJsa-  (20) 

■^0 

where  /?tx  is  the  transmitter  power  (Watts),  gtx  is  the  transmit  antenna  directive  power 
gain  relative  to  isotropic  evaluated  at  the  elevation  angle  of  the  direct  ray  at  the 
transmitter,  4  =\l4iixd  is  the  loss  due  to  signal  spreading,  /ra  is  loss  due  to  rain 
attenuation  and  /ga  is  loss  due  to  gaseous  absorption.  Here  Zo  =  12071  ~  377  Q  is  the 
characteristic  impedance  of  the  atmosphere  and  Xd  is  the  direct  ray  path  length  in  meters 
given  in  equation  (14). 

For  the  reflected  wave  we  first  write 

-:^  =  Pt.gJsKJsa-^  (21) 

■^0 

where  Eq  is  the  field  amplitude  for  a  wave  reflected  from  a  smooth,  flat,  perfectly 
reflecting  surface  as  described  above.  Here  l^\l4itXr  is  the  spreading  loss  for  the 
reflected  wave  path  and  the  other  terms  are  the  same  as  above  with  appropriate  changes 
for  the  different  path.  The  actual  reflected  wave  amplitude  is  obtained  by  inserting  (21) 
into  (18), 


E^ 

-1^  =  PtxgJsKJga^R^D-  ^  (22) 

where  Ir=\R\^  <  1  and  Id=D^  <  1. 

The  total  field  at  the  satellite  is  obtained  from  equations  (16),  (17),  (19),  (21)  and 
(22).  In  these  expressions  the  loss  terms  are  ratios  which  are  less  than  unity  for  a  true 
loss.  To  obtain  4a  from  Zra,  (the  dB  value  derived  in  the  rain  attenuation  routine),  we  use 
/ra=l/10**Tra.  /ga  is  obtained  from  Tga  in  a  similar  manner . 


2.  Rain  Attenuation  Loss 

It  is  well  known  that  absorption  due  to  rain  along  the  propagation  path  is  a  major 
source  of  signal  degradation  for  earth  to  satellite  communications  (book)  and  we  have 
included  a  model  for  estimation  of  rain  loss  in  the  ESPM^  model.  There  exists  a  large 
body  of  research  on  the  estimation  of  the  geographical  and  frequency  dependence  of  rain 
loss  (book)  and  the  International  Telecommunication  Union  (ITU)  periodically  updates 
and  modifies  their  recommended  model  based  on  this  research.  For  the  ESPM^  module 
we  have  chosen  to  implement  the  current  ITU  model  (ITU-R  P.618-8,  P.837-4,  P.838-3, 
P.839-3),  which  was  developed  by  Dissanayake  and  Alnutt  (A.W.  Dissanayake  and  J.E. 
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Alnutt,  1992)  and  has  been  shown  to  provide  accurate  predictions  when  compared  to  data 
collected  over  long  time  frames  over  several  frequencies  (G.S.  Feldhake  and  L.  Ailes- 
Sengers,  2002). 

The  ITU  model  includes  a  contour  map  of  the  maximum  height  above  the  earth  for 
which  rain  is  likely  to  exist  (the  mean  0°  degree  isotherm)  (ITU-R  P.839-3,  2001).  This 
height  is  evaluated  at  the  transmitter  site  and  only  that  portion  of  the  ray  path  lying  below 
the  height  is  included  in  the  loss  determination. 

The  ITU  models  are  statistical  in  nature  and  include  estimates  of  rainfall  rate 
exceeded  for  a  given  probability  of  the  average  year  for  the  transmitter  site  (ITU-R  P- 
837-4,  2003).  For  our  implementation  we  have  chosen  a  probability  level  of  .01,  leading 
to  a  near  ‘worst  case’  estimate  for  long  term  system  planning. 

The  output  of  the  series  of  ITU  models  is  an  estimate  of  the  attenuation  rate  for 
propagation  paths  in  the  region  of  the  transmitter.  We  determine  the  path  loss  by 
multiplying  the  attenuation  rate  by  the  ray  path  length  determined  in  the  ray  optics 
module.  The  path  loss,  for  both  the  direct  and  earth  reflected  rays,  are  then  used  in 
predicting  the  final  signal  strength  as  described  above.  A  detailed  description  of  the  rain 
attenuation  model  is  outside  the  scope  of  this  document  and  the  interested  reader  is 
directed  to  the  ITU  documents  for  more  information. 

3.  Atmospheric  Gaseous  Absorption  Loss 

The  frequency  dependent  loss  of  signal  amplitude  due  to  absorption  by  molecular 
species  in  the  atmosphere  is  another  important  mechanism  that  must  be  included  for 
accurate  prediction  of  signal  strength.  These  losses  are  produced  by  resonant  interactions 
of  the  wave  with  particular  molecular  species  in  the  atmosphere. 

To  calculate  the  gaseous  absorption  for  a  particular  ray,  we  again  turn  to  the  ITU 
recommended  model  (ITU-R  P.676-6,  2005,  Annex  1).  Unlike  the  rain  attenuation 
model,  the  gaseous  absorption  is  a  height  dependent  phenomenon  and  a  gaseous 
attenuation  rate  is  determined  at  each  layer  height  equivalent  to  the  height  at  each  level  of 
refractive  index  input  to  the  ESPM^.  This  data  is  input  to  the  model  along  with  the 
refractivity  profile  and  the  absorption  loss  over  the  entire  propagation  path  is  determined 
for  each  ray  by  summing  the  incremental  losses  within  each  layer.  The  calculation  of 
gaseous  attenuation  rate  is  valid  from  1  to  1000  GHz.  A  detailed  description  of  the 
gaseous  absorption  model  is  outside  the  scope  of  this  document  and  the  interested  reader 
is  directed  to  the  ITU  documents  for  more  information. 

4.  Additional  Output  Parameters:  Phase  Path  Length,  Group  Path  Length  and  Faraday 
Rotation 

Several  parameters  of  general  interest,  in  addition  to  the  signal  strength  and  antenna 
pointing  angles,  can  be  calculated  from  the  ray  trace  results.  These  parameters,  which  are 
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obtained  by  integration  along  the  ray  path,  include  the  phase  path  length,  the  group  path 
length  and  the  rotation  of  the  plane  of  polarization  (Faraday  Rotation). 

In  general,  determination  of  the  signal  strength  at  the  satellite  requires  calculation  of 
the  phase  path  difference  between  the  direct  and  the  earth  reflected  rays  (see  equation 
(19)).  Equation  (17)  defines  the  phase  path  length  'p'  for  a  ray  connecting  transmitter  to 
satellite  as 


P  = 


(23) 


where  ‘f  is  the  transmitter  antenna  location,  ‘s’  is  the  satellite  location  and  ‘rfc’  is  an 
element  of  the  homed  ray  connecting  the  transmitter  and  receiver.  The  total  phase  path  is 
made  up  of  the  sum  of  contributions  from  each  refractive  layer  traversed  by  the  ray  and 
so  we  focus  on  the  contribution  of  a  single  layer  between  r=r,  and  r=r,+/  and  write 


Pi  = 


Vi 


(24) 


Following  the  same  procedure  outlined  earlier,  we  change  the  integral  along  the  ray  path 
to  one  along  the  radial  coordinate  r  (see  equation  (6)).  Then,  making  similar 
approximations  as  in  equations  (7)  and  (8),  we  can  write  the  contribution  to  the  phase 
path  in  the  layer 


Pi  =  nt  +  2n^r^"Y^  +  sf  r^Y 


(25) 


where  Fi,  i=l,2,3,  are  defined  in  equation  (9)  and  the  other  variables  are  also  defined  in 
the  text.  Evaluation  of  Fi,  i=l,2,3,  is  outlined  in  the  APPENDIX.  The  total  phase  path 
length  is  then  obtained  from  the  sum  of  contributions  from  each  layer  traversed.  These 
calculations  are  made  for  both  the  direct  ray  path  and  the  reflected  ray  path  if  the  mode  is 
calculated. 

Another  parameter  which  is  generally  of  interest  for  trans-ionospheric  propagation  is 
the  group  path  delay  or  group  path  length,  which  is  simply  determined  from  the  group 
delay  by  multiplying  by  the  speed  of  light  in  vacuum.  The  group  delay  represents  the 
time  delay  for  the  information  modulated  onto  the  carrier  wave  and  is  especially 
important  for  digital  signals  transmitted  through  the  ionosphere. 

In  equations  (6)  and  (8)  we  defined  the  approximations  used  for  the  atmospheric 
index-of-refraction  in  the  ray  trace  application.  We  note  that  the  index  in  the  atmosphere 
is  not  an  implicit  function  of  the  frequency  of  the  wave.  In  that  case,  the  group  path 
length  and  the  phase  path  length  are  identical.  However,  an  electromagnetic  wave  that 


22 


traverses  any  part  of  the  earth’s  ionosphere  experiences  additional  delay  and  the  phase 
path  length  and  group  path  length  will  differ.  To  derive  an  expression  for  group  delay 

VI.  APPENDIX 


In  the  text  we  claimed  that  equation  (5),  the  basic  ray  equation,  can  be  written  with 
equations  (7)  and  (8)  as 


2  3 

r,—r,+-r, 

K  r, 


(Al) 


The  definitions  of  the  Ti  are  given  in  equation  (9)  and  the  form  of  equation  (Al)  is 
obvious  from  equation  (7).  In  this  appendix  we  wish  to  outline  the  derivation  of  the 
solution  of  the  ray  equation  given  in  equation  (10). 

We  first  note  that  the  integrals  r2  and  can  be  written  in  terms  of  Ti  as  follows  (R.S. 
Burington,  1965) 


c  2c 

r,  (h)  =  ^  D{h)  -  ^  r,  (A)  -  ^  r,  (A) 
2c  4c  2c 


where  D(h)={a+bh+ch^y^  and  a,b,c  are  defined  in  the  text  after  equation  (9).  Inserting 

(A3) 


(A2) 


r2  and  Ta  into  equation  (Al),  and  with  equation  (A2),  we  write  (Al)  as 

^/.|:;-  =  Cj^r,(A)  +  Z)(A)5(A)]5» 


where  A  and  B(h)  are  defined  in  the  text  after  equation  (10).  For  the  case  c  <  0,  the 
integral  T i  is  given  by 


-4ab 


and  for  c  >  0  we  have 

rj(A)  =  -plog(2cA  +  A  +  l^fc-^a  +  bh  +  ch^ ) . 
ylc 

Note  that  T i(A)=F(A)  in  the  text.  Evaluation  of  equation  (A3)  at  the  integration  boundaries 
results  in  equation  (10)  in  the  text. 
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For  the  special  case  c  =  0,  the  integrals  ¥2  and  Fs  can  again  be  written  in  terms  of  Fi, 
which  is  given  in  this  case  by  F \{h)={2lb){a+bh)^'^ .  Evaluation  of  equation  (A2)  leads  to 


or 


v\T-c„\Q(h)r,{h)l 


(A4) 


where  Q{h)=\+2{2a-bh)/3bri+{Sa^-Aabh+3b^h^)l5b^r^ . 

Finally,  we  note  that  for  a  down-going  ray  incident  on  the  refractive  layer  (^i+i  >  90°) 
the  form  of  these  equations  remains  the  same.  However,  the  signs  of  the  r2  and  T3  terms 
in  equation  (Al)  are  reversed.  Also,  for  the  coefficients  b  and  c  (see  text  after  equation 
(9)),  the  sign  of  the  second  term  in  each  is  reversed.  Further  we  note,  in  this  case,  c  is 
always  positive. 
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